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1  SUMMARY 


Supported  by  AFOSR,  our  previous  direct  numerical  simulation  (DNS)  studies  from  December 
1999  to  November  2002  have  demonstrated  that  our  unique  high-order  shock  fitting  simulation 
approach  is  a  powerful  tool  in  studying  supersonic  and  hypersonic  boundary-layer  stability  and 
transition  physics.  However,  experimental  validation  of  numerical  simulation  methods  for 
hypersonic  boundary  layer  stability  and  transition  has  not  been  done  and  is  urgently  needed.  The 
objective  of  this  research  project  is  to  compare  our  numerical  simulation  solutions  with  available 
experimental  or  theoretical  results  on  hypersonic  boundary  layer  receptivity  and  stability;  and  to 
conduct  extensive  DNS  studies  on  the  flow  mechanisms  of  hypersonic  boundary  layer  receptivity 
and  stability.  During  the  three-year  period,  we  have  conducted  extensive  DNS  studies  on  the 
receptivity  of  hypersonic  boundary  layer  flows  over  a  sharp  wedge,  a  flat  plate,  a  blunt  cone,  and 
the  FRESH  aeroshell.  DNS  studies  are  compared  with  Stetson's  1984  stability  experiment  on  Mach 
7.99  flow  over  a  blunt  cone,  and  Maslov’s  leading-edge  receptivity  experiment  on  Mach  5.92  flow 
over  a  flat  plate.  We  also  collaborate  with  Prof.  Tumin  in  University  of  Arizona  to  compare 
numerical  and  theoretical  results  on  receptivity  of  a  Mach  8  flow  over  a  sharp  wedge  to  wall 
blowing-suction,  and  to  deeply  analyze  the  receptivity  characteristics.  Our  numerical  studies  have 
been  validated  to  be  of  high  accuracy  and  led  to  further  understanding  of  hypersonic  boundary  layer 
receptivity  mechanism.  Such  understanding  can  lead  to  better  tool  for  the  prediction  and  control  of 
high  speed  boundary  layer  transition. 

The  main  research  contributions  are: 

1.  We  have  conducted  DNS  studies  on  the  receptivity  of  a  Mach  8.0  flow  over  a  sharp  wedge  of  a 
half-angle  of  5.3°  to  wall  blowing-suction.  The  results  show  that  Mode  F,  mode  S,  and  acoustic 
modes  are  simultaneously  excited  by  wall  blowing-suction  disturbances.  It  is  found  that  the 
frequency  of  the  blowing-suction  actuator  has  a  significant  effect  on  the  receptivity  process 
owing  to  the  frequency  dependence  of  the  synchronization  point  between  mode  F  and  mode  S. 
Furthermore,  the  excitation  of  mode  S  is  strongly  affected  by  profile  and  length  of  the  forcing 
actuator.  All  cases  of  numerical  simulations  consistently  show  that  the  synchronization  point 
plays  an  important  role  in  the  excitation  of  mode  S  by  wall  blowing-suction,  i.e.,  mode  S  is 
strongly  excited  only  when  the  blowing-suction  actuator  is  located  upstream  of  the 
synchronization  point.  When  the  forcing  actuator  is  downstream  of  the  synchronization  point, 
there  is  very  little  excitation  of  mode  S,  despite  the  fact  that  the  blowing-suction  actuator  is  still 
located  within  the  unstable  region  of  mode  S. 

2.  We  collaborated  with  Prof.  Tumin  in  University  of  Arizona  to  compare  numerical  and 
theoretical  results  on  receptivity  of  the  Mach  8  flow  over  a  sharp  wedge  to  wall  blowing-suction, 
and  to  deeply  analyze  the  receptivity  characteristics.  The  perturbation  flow  field  downstream  of 
the  blowing-suction  actuator  was  decomposed  into  boundary-layer  wave  modes  with  the  help  of 
the  biorthogonal  eigenfunction  system.  It  was  shown  that  amplitudes  of  boundary-layer  wave 
modes  calculated  with  the  help  of  the  theoretical  receptivity  model  had  good  agreements  with 
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those  obtained  by  projecting  the  numerical  results  onto  the  boundary-layer  wave  modes. 

3.  We  have  studied  the  receptivity  of  a  Mach  5.92  over  a  flat  plate  to  various  two-dimensional  wall 
perturbations:  oscillation,  blowing-suction,  and  temperature  perturbation.  The  numerical  results 
show  that  all  types  of  two-dimensional  wall  perturbations  eventually  result  in  the  same  type  of 
instability  wave  (mode  S)  in  the  boundary  layer.  The  hypersonic  boundary-layer  flow  is  most 
sensitive  to  wall  blowing-suction  and  least  sensitive  to  wall  temperature  perturbation.  It  is 
noticed  that  the  profile  of  temperature  perturbation  has  a  strong  effect  on  receptivity  process. 
The  receptivity  of  the  hypersonic  boundary  layer  to  temperature  perturbation  is  proportional  to 
the  energy  introduced  to  the  steady  base  flow,  where  the  energy  is  determined  by  the  profile  of 
temperature  perturbation. 

4.  We  have  conducted  DNS  studies  on  the  transient  growth  phenomena  of  a  Mach  5.92  flat-plate 
boundary  layer  and  a  Mach  7.99  flow  over  a  blunt  cone  to  small-scale  three-dimensional 
stationary  roughness  element.  The  surface  roughness  is  periodic  in  spanwise  direction.  Effect  of 
thermal  flat  plate  boundary  conditions  on  the  receptivity  process  is  considered  by  comparing 
numerical  simulations  on  adiabatic  and  isothermal  flat  plates,  respectively.  The  numerical 
results  show  that  counter  rotating  streamwise  vortices  are  excited  by  stationary  roughness 
element.  The  temperature  increases  in  the  wake  region  of  roughness  peak  whereas  it  decreases 
in  the  wake  region  of  roughness  valley.  Furthermore,  the  speed  of  wake  from  roughness  peak  is 
slower  than  that  from  roughness  valley.  All  these  results  are  consistent  with  those  of  Tumin's 
theoretical  analysis.  It  is  also  concluded  based  on  the  numerical  results  that  the  roughness 
element  on  adiabatic  flat  plate  is  more  efficient  in  the  excitations  of  streamwise  vortices  and 
transient  growth. 

5.  We  developed  a  new  high-order  immersed  interface  method  for  computing  boundary  layer 
transition  with  surface  roughness.  The  work  was  motivated  by  discussions  among  the 
participants  in  the  Transition  Study  Group  Open  Forum  Held  in  the  AIAA  meeting  in  Reno  in 
January  2005.  It  was  pointed  out  that  that  there  is  a  need  to  compute  boundary  layer  stability 
and  transient  growth  with  distributive  surface  roughness.  The  new  method  can  be  arbitrarily 
high-order  accuracy  in  the  whole  flow  field,  including  the  interface  with  discontinuity,  and  it 
can  be  useful  for  the  simulations  of  transient  growth  to  finite  or  nonlinear  surface  roughness. 


2  RESEARCH  OBJECTIVES 

The  prediction  of  laminar-turbulent  transition  in  supersonic  and  hypersonic  boundary  layers  is  a 
critical  part  of  the  aerodynamics  and  heating  analyses  for  the  development  of  hypersonic 
transportation  vehicles  and  re-entry  vehicles.  The  key  for  an  accurate  transition  prediction  is  the 
understanding  of  the  physical  mechanisms  that  lead  to  transition.  However,  many  important 
physical  mechanisms  leading  to  hypersonic  boundary  layer  transition  are  currently  still  not  well 
understood.  Among  them,  receptivity,  which  can  provide  the  missing  link  between  environment 
disturbances  and  boundary  layer  instability  growth  and  transition,  is  currently  poorly  understood. 
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Due  to  the  difficulties  in  conducting  hypervelocity  experiments  and  the  complexity  of  hypersonic 
flows,  fundamental  hypersonic  studies  increasingly  rely  on  the  use  of  DNS  as  a  research  tool.  In 
recent  years,  DNS  has  become  a  powerful  tool  in  the  studies  of  the  stability  and  receptivity  of 
supersonic  and  hypersonic  boundary  layers. 

Supported  by  AFOSR,  we  have  developed  high-order  shock-fitting  DNS  methods  [1],  which  can  be 
directly  applied  to  hypersonic  boundary  layers  over  realistic  blunt  bodies  with  the  effects  of  nose 
bluntness,  the  presence  of  bow  shock  waves,  and  the  real-gas  effects  at  high  temperatures.  In  the 
past  several  years,  we  have  studied  the  receptivity  and  stability  of  a  number  of  2-D  and  3-D 
hypersonic  flows  over  sharp  wedges  [2,  3],  .flat  plates  [4,  5],  and  blunt  cones  [6-8]. 

The  main  objectives  of  this  research  are  to  compare  our  numerical  simulation  solutions  with 
available  experimental  or  theoretical  results  on  hypersonic  boundary  layer  receptivity  and  stability, 
and  to  conduct  extensive  DNS  studies  on  the  flow  mechanisms  of  hypersonic  boundary  layer 
receptivity  and  stability.  The  DNS  numerical  tools,  as  well  as  other  theoretical  approaches  such  as 
the  linear  stability  analysis  and  the  modal  decomposition  of  biorthogonal  eigenfunction  system,  are 
used  to  gain  new  fundamental  understandings  of  receptivity  and  stability  of  hypersonic  boundary 
layers.  The  goal  of  these  studies  is  to  lead  to  a  better  understanding  of  validity,  accuracy,  and 
limitation  of  DNS  results  through  experimental  and  theoretical  comparisons,  and  to  achieve  new 
improvements  on  the  prediction  tools  for  hypersonic  boundary  layer  transition,  by  incorporating  the 
effects  of  forcing  disturbances. 

Six  major  research  tasks  have  been  undertaken  in  this  research  project.  They  are: 

A)  Code  validation  of  Stetson’s  stability  experiment  of  Mach  7.99  flow  over  a  blunt  cone  [9]; 

B)  Receptivity  to  wall  blowing-suction  of  a  Mach  8  flow  over  a  sharp  wedge; 

C)  Receptivity  to  various  two-dimensional  wall  perturbations  of  the  flat-plate  boundary  layer  flow 
corresponding  to  Maslov  et  al.’s  leading-edge  receptivity  experiments  [10]; 

D)  Roughness  induced  transient  growth  in  a  Mach  7.99  flow  over  a  blunt  cone  and  the  Mach  5.92 
flat-plate  boundary  layer; 

E)  Development  of  a  new  high-order  immersed  interface  method  [11,  12]  for  the  simulations  of 
transient  growth  to  finite  or  nonlinear  surface  roughness; 

F)  Stability  and  receptivity  of  hypersonic  flow  over  the  FRESH  areoshell. 

The  main  approach  is  to  use  DNS  as  a  research  tool  to  study  hypersonic  boundary  layer  receptivity 
physics.  Linear  stability  analysis  (LST)  and  model  decomposition  of  biorthogonal  eigenfunction 
system  are  used  in  conjunction  with  DNS  to  help  better  understand  the  flow  instability  mechanisms 
in  numerical  solutions.  Our  fifth-order  upwind  finite  difference  shock  fitting  method  [13]  for  the 
DNS  of  hypersonic  flows  with  a  strong  bow  shock  is  used.  The  use  of  the  high-order  shock-fitting 
scheme  makes  it  possible  to  obtain  highly  accurate  mean  flow  and  unsteady  solutions,  which  are 
free  of  spurious  numerical  oscillations  behind  the  bow  shock. 


3  RESEARCH  ACCOMPLISHMENTS 
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The  research  project  described  in  the  previous  section  was  carried  out  in  a  three  period  from 
January  1,  2004  to  January  1,  2007.  The  three-year  research  project  supported  by  this  grant  has  led 
to  the  publications  of 

1 .  Nine  papers  in  the  archive  journals  or  book  chapters, 

2.  Thirteen  conference  papers  (mainly  AIAA  papers), 

3.  One  Ph.D.  thesis  by  a  student  who  has  graduated  with  his  Ph.D.  degree. 

More  importantly,  the  research  has  produced  a  large  amount  of  new  results,  and  has  led  to  new 
understanding  of  a  number  of  hypersonic  boundary  layer  receptivity  mechanisms.  The 
understanding  of  the  receptivity  process  is  the  key  to  improve  current  hypersonic  boundary  layer 
transition  prediction  methods.  In  the  following  sections,  the  major  research  accomplishments  are 
summarized,  followed  by  a  more  details  discussions  of  some  of  our  new  results  in  hypersonic 
boundary  layer  receptivity. 

3.1  Code  validation  of  Stetson’s  stability  experiment  of  Mach  7.99  flow  over  a  blunt  cone 
Publications:  [4],  [15],  [18] 

We  have  carried  out  code  validation  for  the  Stetson’s  stability  experiment  of  Mach  7.99  flow  over  a 
blunt  cone.  The  simulation  results  for  the  second  mode  growth  compared  well  with  those  of  the 
experiments  and  LST.  The  receptivity  to  wall  blowing/suction  at  three  different  surface  locations  is 
used  to  excite  the  second  mode  instability  waves  in  the  numerical  simulation  of  the  full 
Navier-Stokes  equations.  It  is  found  (Fig.  1)  that  the  experimental  second-mode  amplitude  growth 
ratios,  N  factors,  and  growth  rates  agree  well  with  those  of  the  simulation  in  a  linear  growth  region 
of  195<s<215  for  frequencies  near  the  peak  second  mode  instability.  It  is  concluded  that  the 
simulation  and  LST  are  accurate  in  predicting  second  mode  growth  in  the  current  hypersonic 
boundary  layer  flows. 


Fig.  1  The  7°  half-angle  blunt  cone  of  Stetson's  Mach  7.99  stability  experiments  [9]  and  numerical  simulation 
of  second  mode  growth  by  surface  blowing/suction  at  surface  station  of  215  nose  radii. 
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Fig.  2  Experimental,  simulation  and  LST  comparison  of  growth  rates,  wave  numbers,  and  amplitude  ratios 
(using  195  surface  location  as  reference)  at  surface  station  of  s  =  215. 


Fig.  2  compares  the  second  mode  growth  rates,  and  wave  numbers  at  the  surface  station  of  5=215. 
Two  sets  of  experimental  growth  rates  of  Stetson  are  plotted.  The  first  set  (labeled  as  Stetson's 
Experiment  II  in  the  figure)  is  the  growth  rates  published  in  Stetson’s  original  1984  paper.  The 
second  set  (Stetson's  Experiment  I)  are  those  calculated  by  the  current  author  based  on  the  new 
published  experimental  results  (Fig.  28  of  Schneider's  review  paper  [14]).  The  growth  rates  of  the 
current  simulation  of  four  frequencies  are  plotted  in  the  figure.  These  four  frequencies  show  second 
mode  excitation  at  this  location  in  the  simulation  results.  This  figure  shows  that  the  growth  rates  of 
the  full  Navier-Stokes  simulation  agree  very  well  with  the  LST  results.  On  the  other  hand,  the  newly 
published  experimental  growth  rates  are  much  closer  to  those  of  the  original  LST  and  experimental 
results  with  a  3%  difference  in  wave  number  and  a  12%  difference  in  growth  rate  for  the  peak 
growth  rate.  This  is  a  reasonable  agreement  considering  the  fact  that  the  wall  surface  is  not  perfectly 
adiabatic  in  the  experiment. 

3.2  Receptivity  to  wall  blowing-suction  of  a  Mach  8  flow  over  a  sharp  wedge 
Publications:  [3],  [7],  [14],  [16],  [21] 

Understanding  the  receptivity  process  is  of  critical  importance  to  the  prediction  and  control  of 
boundary  layer  transition.  With  the  development  of  advanced  computers  and  numerical  techniques, 
numerical  simulation  of  the  receptivity  process  by  directly  solving  Navier-Stokes  equations  has 
become  feasible.  By  solving  the  compressible  linearized  Navier-Stokes  equations,  Malik  et  al.  [15] 
investigated  the  responses  of  a  Mach  8  flow  over  a  sharp  wedge  of  a  half-angle  of  5.3°  to  three 
types  of  external  forcing:  a  planar  free-stream  acoustic  wave,  a  narrow  acoustic  beam  enforced  on 
the  bow  shock  near  the  leading  edge,  and  a  blowing-suction  slot  on  the  wedge  surface.  They 
concluded  that  these  three  types  of  forcing  eventually  resulted  in  the  same  type  of  instability  waves 
in  the  boundary  layer.  However  the  receptivity  mechanism  was  not  studied  in  detail.  Ma  and  Zhong 
[16]  studied  the  receptivity  mechanisms  of  the  same  hypersonic  boundary  layer  to  various 
freestream  disturbances,  i.e.,  fast  and  slow  acoustic  waves,  vorticity  waves,  and  entropy  waves,  by 
solving  the  two-dimensional  compressible  Navier-Stokes  equations.  They  found  that  the  stable 
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modes  in  the  boundary  layer  played  a  very  important  role  in  the  receptivity  process. 

In  this  project,  the  receptivity  mechanism  of  the  Mach  8  flow  over  the  sharp  wedge  to  wall 
blowing-suction  is  investigated  in  detail  by  series  of  numerical  simulations.  The  steady  base  flow  is 
simulated  by  solving  the  two-dimensional  compressible  Navier-Stokes  equations  with  a 
combination  of  a  fifth-order  shock-fitting  finite  difference  method  and  a  second-order  TVD  scheme. 
In  the  receptivity  simulations,  periodic  blowing-suction  disturbances  are  introduced  to  the  steady 
base  flow  through  a  forcing  actuator  on  the  wedge  surface.  The  characteristics  of  boundary-layer 
wave  modes  are  identified  and  evaluated  by  comparing  the  results  of  LST  and  numerical 
simulations.  The  effects  of  frequency  and  location  of  the  forcing  actuator  on  the  receptivity  process 
are  studied  by  considering  seven  cases  of  different  blowing-suction  locations  with  15  frequencies 
for  each  case.  Furthermore,  the  effects  of  profile  and  length  of  the  forcing  actuator  are  investigated. 
For  all  cases  considered  in  this  paper,  forcing  waves  are  weak  enough  so  that  the  disturbances  in  the 
boundary  layer  are  linear.  The  contours  of  the  instantaneous  pressure  perturbation  induced  by 
blowing-suction  disturbance  are  shown  in  Fig.  3.  After  the  blowing-suction  slot,  the  excited 
pressure  perturbations  are  divided  into  two  branches:  one  branch  penetrates  the  boundary  layer  and 
propagates  along  the  Mach  lines  (acoustic  waves)  while  the  other  branch  stays  within  the  boundary 
layer,  which  is  the  unstable  mode  (mode  S),  indicated  by  the  typical  wave  structures  near  the  wall. 
The  amplification  of  instantaneous  pressure  perturbation  along  the  wedge  surface  also  indicates  the 
excitation  of  mode  S. 


Fig.  3  Numerical  simulation  of  receptivity  of  the  Mach  8  flow  over  a  sharp  wedge  to  wall  blowing-suction:  1) 
instantaneous  pressure  perturbation  field;  2)  instantaneous  pressure  perturbation  along  the  surface. 

The  numerical  results  show  that  mode  F,  mode  S,  and  acoustic  modes  are  simultaneously  excited  by 
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wall  blowing-suction  disturbances.  Downstream  of  the  forcing  region,  mode  F  decays  owing  to  its 
inherent  stability  whereas  mode  S  grows  because  of  its  instability.  Acoustic  modes  radiate  into  the 
external  flow  outside  the  boundary  layer.  As  a  result,  mode  S  eventually  becomes  the  dominant 
wave  mode  in  the  boundary  layer.  It  is  found  that  the  frequency  of  the  blowing-suction  actuator  has 
a  significant  effect  on  the  receptivity  process  owing  to  the  frequency  dependence  of  the 
synchronization  point  between  mode  F  and  mode  S.  The  numerical  results  also  show  that  the 
excitation  of  mode  S  is  strongly  affected  by  profile  and  length  of  the  forcing  actuator.  All  cases  of 
numerical  simulations  consistently  show  that  the  synchronization  point  plays  an  important  role  in 
the  excitation  of  mode  S  by  wall  blowing-suction,  i.e.,  mode  S  is  strongly  excited  only  when  the 
blowing-suction  actuator  is  located  upstream  of  the  synchronization  point.  When  the  forcing 
actuator  is  downstream  of  the  synchronization  point,  there  is  very  little  excitation  of  mode  S,  despite 
the  fact  that  the  blowing-suction  actuator  is  still  located  within  the  unstable  region  of  mode  S. 

a)  b) 
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Fig.  4  Comparison  of  the  theoretical  prediction  for  the  receptivity  coefficient  with  data  filtered  out  from  the 
computational  results  for  receptivity  of  a  hypersonic  boundary  layer  to  wall  blowing-suction:  a)  mode  S  and  b) 

mode  F. 

On  the  other  hand,  numerical  simulation  results  require  careful  analytical  analysis  in  order  to  gain  a 
deeper  understanding  of  the  underlining  physics  in  high-speed  boundary-layer  transition.  A  close 
collaboration  between  numerical  and  theoretical  study  of  the  receptivity  phenomena  will  not  only 
verify  both  sets  of  results,  but  also  lead  to  better  understanding  of  flow  physics.  We  have  started 
working  in  this  direction  with  Professor  Tumin  [2,  17]  in  University  of  Arizona.  Specifically,  direct 
numerical  simulation  of  receptivity  in  a  boundary  layer  over  a  sharp  wedge  of  half-angle  5.3°  was 
carried  out  with  two-dimensional  perturbations  introduced  into  the  flow  by  periodic-in-time 
blowing-suction  through  a  slot.  The  free-stream  Mach  number  was  equal  to  8.  The  perturbation 
flow-field  downstream  from  the  slot  was  decomposed  into  normal  modes  with  the  help  of  the 
biorthogonal  eigenfunction  system.  Filtered-out  amplitudes  of  two  discrete  normal  modes  and  of  the 
fast  acoustic  modes  are  compared  with  the  linear  receptivity  problem  solution.  Fig.  4  shows  a 
comparison  of  the  theoretical  receptivity  coefficient  with  the  amplitude  filtered  out  from  the 
computational  results  by  the  theoretical  analysis.  One  can  see  from  the  figure  that  there  is  an 
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excellent  agreement  between  amplitudes  calculated  with  the  help  of  the  theoretical  receptivity 
model. 

3.3  Receptivity  to  various  two-dimensional  wall  perturbations  of  a  Mach  5.92  flat-plate 
boundary  layer  flow 

Publications:  [10],  [13] 

The  receptivity  of  a  Mach  5.92  boundary-layer  flow  over  a  flat  plate  to  periodic  two-dimensional 
wall  perturbations  is  investigated  by  numerical  simulations  and  linear  stability  theory  (LST).  The 
flow  conditions  are  the  same  as  those  of  Maslov  et  al.’s  leading-edge  experiment  [10].  The  steady 
base  flow  is  simulated  by  solving  two-dimensional  compressible  Navier-Stokes  equations  with  a 
combination  of  a  fifth-order  shock-fitting  method  and  a  second-order  TVD  scheme.  Accuracy  of  the 
numerical  steady  base  flow  is  validated  by  comparing  with  the  theoretical  self-similar 
boundary-layer  solution  and  the  published  experimental  results.  The  characteristics  of 
boundary-layer  wave  modes  are  identified  and  evaluated  by  comparing  the  results  of  LST  and 
numerical  simulations.  In  receptivity  simulations,  three  types  of  periodic  two-dimensional  wall 
perturbations,  oscillation,  blowing-suction,  and  temperature  perturbation,  are  considered.  All  wall 
perturbations  are  introduced  to  the  steady  base  flow  by  a  forcing  slot  located  on  the  flat  plate. 
Effects  of  the  perturbation  type  and  the  profile  of  temperature  perturbation  on  receptivity  process 
are  studied  by  considering  six  cases  of  wall  perturbations.  Fig.  5  compare  normalized  Mach  number 
and  dimensionless  streamwise  velocity  distributions  across  the  boundary  layer  at  three  different 
locations  of  x*  =  96  mm,  121  mm,  and  138  mm.  The  numerical  results  agree  well  with  the 
experimental  results  and  the  boundary-layer  solution  near  the  plate.  However,  in  the  region  of  n  > 
5,  the  numerical  results  have  a  better  agreement  with  the  experimental  results.  The  difference 
between  the  numerical  results  and  the  boundary-layer  solution  is  mainly  caused  by  the  existence  of 
the  bow  shock,  because  the  effect  of  the  bow  shock  is  neglected  in  the  calculation  of  the 
compressible  boundary-layer  equations. 


Fig.  5  Numerical  simulation  of  a  Mach  5.92  flat-plate  boundary  layer  to  two-dimensional  wall  perturbations  and 
comparisons  of  normalized  Mach  number  and  dimensionless  streamwise  velocity. 
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Fig.  6  Numerical  simulation  of  receptivity  of  the  Mach  5.92  flow  over  a  flat  plate  to  two-dimensional  wall 
perturbations:  1)  instantaneous  pressure  perturbation  along  the  wall;  2)  comparison  of  numerical  and  LST 
growth  rates;  3)  comparison  of  numerical  and  LST  wave  number;  4)  amplitude  comparison  of  pressure 
perturbation  (100  kHz)  for  cases  of  wall  oscillation,  blowing-suction,  and  temperature  perturbation;  5) 
amplitude  comparison  of  pressure  perturbation  (100  kHz)  for  temperature  perturbations  of  different  profile. 

Fig.  6  shows  that  all  six  cases  of  two-dimensional  wall  perturbations  eventually  result  in  the  same 
type  of  instability  wave  (mode  S)  in  the  boundary  layer.  The  hypersonic  boundary-layer  flow  is 
most  sensitive  to  wall  blowing-suction  and  least  sensitive  to  wall  temperature  perturbation.  It  is 
noticed  that  the  profile  of  temperature  perturbation  has  a  strong  effect  on  receptivity  process.  The 
receptivity  of  the  hypersonic  boundary  layer  to  temperature  perturbation  is  proportional  to  the 
energy  introduced  to  the  steady  base  flow,  where  the  energy  is  determined  by  the  profile  of 
temperature  perturbation. 

3.4  Roughness  induced  transient  growth  in  a  Mach  7.99  flow  over  a  blunt  cone  and  the  Mach 
5.92  flat-plate  boundary  layer 

Publications:  [10],  [11],  [15] 

Recent  theoretical  studies  have  shown  the  importance  of  transient  growth  and  surface  roughness  in 
bypass  transition  [18,  19].  Reshotko  and  Tumin  [20]  have  developed  a  model  for  roughness  induced 
transient-growth  transition.  So  far,  there  has  not  been  any  direct  numerical  simulation  study  on 
transient  growth  in  hypersonic  boundary  layers.  Furthermore,  it  is  not  known  how  the  optimal 
disturbances  computed  by  the  transient  growth  theory  are  generated  by  surface  roughness.  These  are 
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important  issues  related  to  the  receptivity  of  transient  growth.  We  use  direct  numerical  simulation  to 
conduct  such  studies  because  the  complete  receptivity  process  and  transient  growth  can  be 
simulated  with  minimum  simplification  assumptions.  In  the  review  period,  we  have  done  simulation 
studies  on  roughness  induced  transient  growth  for  hypersonic  flow  over  blunt  bodies  and  a  flat  plate. 
The  basic  approach  is  to  impose  stationary  3-D  surface  roughness  in  the  early  stable  region  of  the 
steady  base  hypersonic  boundary  layers.  The  subsequent  development  of  disturbances  in  the 
boundary  layer  is  captured  by  numerical  simulations. 


Fig.  7  Simulation  solutions  of  temperature  perturbations  due  to  distributed  surface  roughness  for  Case  D  of 
Mach  7.99  flow  over  a  blunt  cone  with  88.9  mm  nose  radius. 


The  specific  flow  conditions  for  the  blunt  cones  are  the  same  as  Stetson’s  Mach  7.99  flow  over  a  7 
deg.  half  angle  blunt  cones  with  the  four  sets  of  nose  radii  (Cases  A  to  D).  So  far,  only  small 
distributed  surface  roughness  has  been  imposed  on  the  surface  in  the  following 

form:  h(s,  0)  =  sh0  An  cos(2 nnO  +  (j)  .  Fig.  7  shows  the  results  of  temperature  perturbations 

n 

for  Case  D  with  nose  radius  of  88.9  mm.  The  flow  perturbation  induced  by  weak  surface  roughness 
shows  the  growth  of  disturbances  along  the  cone  surface.  Work  is  currently  under  way  to  extend 
these  initial  studies  to  longer  computational  domain  and  to  study  the  transient  growth  properties 
induced  by  surface  roughness. 


Fig.  8.  Numerical  simulation  of  receptivity  of  the  Mach  5.92  flow  over  a  flat  plate  to  three-dimensional 
stationary  roughness  element:  1)  vector  plot  of  wall-normal  and  spanwise  velocity  perturbations  in  the  wake 
region;  2)  contours  of  temperature  perturbation  in  the  wake  region;  3)  contours  of  streamwise  velocity 

perturbation  in  the  wake  region. 

We  have  also  studied  the  receptivity  to  three-dimensional  stationary  roughness  elements  of  a 
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hypersonic  flat-plate  boundary  layer,  corresponding  to  Maslov  et  al.’s  leading-edge  experiment  [10]. 
The  surface  roughness  is  periodic  in  spanwise  direction.  Effect  of  thermal  flat  plate  boundary 
conditions  on  the  receptivity  process  is  considered  by  comparing  numerical  simulations  on  adiabatic 
and  isothermal  flat  plates,  respectively.  Fig.  8  shows  that  counter  rotating  streamwise  vortices  are 
excited  by  stationary  roughness  element.  The  temperature  increases  in  the  wake  region  of  roughness 
peak  whereas  it  decreases  in  the  wake  region  of  roughness  valley.  Furthermore,  the  speed  of  wake 
from  roughness  peak  is  slower  than  that  from  roughness  valley.  It  is  also  concluded  based  on  the 
numerical  results  that  the  roughness  element  on  adiabatic  flat  plate  is  more  efficient  in  the 
excitations  of  streamwise  vortices  and  transient  growth. 

3.5  Development  of  a  new  high-order  immersed  interface  method 

Publications:  [2],  [12] 

During  three-year  project,  we  developed  a  new  high-order  immersed  interface  method  for 
computing  boundary  layer  transition  with  surface  roughness.  The  work  was  motivated  by 
discussions  among  the  participants  in  the  Transition  Study  Group  Open  Forum  Held  in  the  AIAA 
meeting  in  Reno  in  January  2005.  It  was  pointed  out  that  that  there  is  a  need  to  compute  boundary 
layer  stability  and  transient  growth  with  distributive  surface  roughness.  The  methods  for  simulating 
flow  with  irregular  surface  are  similar  to  those  for  simulation  of  multi-phase  flows  or  fluid-structure 
interaction.  A  major  difficulty  in  such  simulation  is  due  to  the  fact  that  the  flow  variables  and  their 
derivatives  are  not  continuous  across  the  interface  or  irregular  surface.  One  way  to  simulate  surface 
roughness  of  an  arbitrary  shape  is  to  use  the  immersed  boundary  method  [21]  or  the  immersed 
interface  method  (IIM)  of  Laveque  and  Li  [22]  with  a  fixed,  non-body  fitted  grid.  The  main 
drawback  of  the  immersed  boundary  method  is  the  fact  that  it  is  at  most  first  order  accurate  at  the 
interface.  The  immersed  interface  method  can  achieve  a  second  order  accuracy.  But  the  IIM  method 
is  difficult  to  be  used  in  complex  flow  problems.  In  order  to  maintain  a  second  order  accuracy,  it  is 
necessary  to  know  the  jump  conditions  of  flow  variables  and  the  jump  conditions  for  the  first  and 
the  second  derivatives  at  the  interface. 


For  a  uniformly  high-order  simulation  of  flow  with  surface  roughness,  we  have  developed  a  new 
high-order  IIM  method  [11].  The  new  method  can  be  arbitrarily  high-order  accuracy  in  the  whole 
flow  field,  including  the  interface  with  discontinuity.  So  far,  the  new  uniformly  2nd  and  4th  order 
accurate  IIM  methods  have  been  tested  for  a  number  of  1-D  and  2-D  model  equations  with 
discontinuous  interfaces.  Figure  5  shows  the  result  of  a  2-D  example  used  by  Leveque  and  Li  [22] 


for  the  following  elliptical  equation, 


d2u  d2u 
dx2  dy2 


^2S{x- X(s))5(y  -Y(s))ds ,  where  the  interface 


is  a  circle  of  radius  0.5  centered  at  the  origin.  There  is  a  jump  for  the  normal  velocity  gradient  at  the 
interface.  A  uniform  grid  of  N  points  is  used  in  the  numerical  computations.  The  results  are 
compared  for  three  different  numerical  methods:  (1)  our  IIM  method  A:  2nd  order  accuracy  in  both 
regular  and  interface  grids,  (2)  our  IIM  method  B:  4th  order  accuracy  in  regular  grids  and  3rd  order  in 
irregular  grid  points,  and  (3)  Leveque  and  Li’s  2nd  order  IIM  [22],  Fig.  9  shows  that  the  new  IIM 
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method  captures  the  interface  accurately.  Table  1  shows  that  the  new  4th  order  IIM  method  is  most 
accurate  among  these  three  methods.  The  results  of  this  test  case  and  other  test  cases  show  the 
potential  of  the  new  high-order  IIM  method  in  accurately  computing  flow  with  irregular  surface,  as 
well  as  two-phase  flow  problems.  We  plan  to  apply  the  new  methods  to  the  simulation  of  boundary 
layer  transient  growth  with  distributed  surface  roughness. 


Table  1 .  Comparison  of  Numerical  Solutions  of  Three  Methods 


Leveque  and  Li’s  IIM 
method  (2nd  order)  [4] 

Current  IIM  method  A 
(2nd  order) 

Current  IIM  method  B 
(4th  order) 

Errors  (80  by  80  grid) 

2.45  xKT4 

6.70  xlO-5 

2.31xl0"6 

Error  Ratios 

3.41 

3.90 

12.6 

Fig.  9  Solution  contours  with  discontinuous  gradient  at  the  interface  and  comparison  of  exact  solution  with  the 
new  2nd  order  IIM  method  A  for  computing  the  2-D  model  equation  with  a  singular  source  term. 


3.6  Stability  and  receptivity  of  hypersonic  flow  over  the  FRESH  areoshell 

At  the  request  of  Dr.  John  Schmisseur,  we  worked  with  Drs.  Datta  Gaitonde  and  Roger  Kimmel  of 
AFRL  on  a  computational  study  of  the  stability  and  receptivity  of  the  FRESH  aeroshell.  One  of  the 
objectives  of  the  FRESH  is  to  measure  smooth-body  and  roughness-induced  transition  on  the  cone 
surface.  The  aeroshell  has  a  spherical  cone  with  a  7  degree  half  angle.  It  uses  either  a  sharp  nose  or 
a  slightly  blunt  nose.  The  objectives  of  our  DNS  studies  are  to  compare  our  results  with  other 
simulation  results  obtained  by  AFRL  and  to  do  DNS  studies  of  the  stability  and  receptivity  of 
hypersonic  boundary  layers  over  the  aeroshell.  So  far,  we  have  finished  the  computations  of  the 
mean  flow  of  a  test  case  with  the  same  conditions  as  a  wind  tunnel  test  of  the  model.  The  flow 

conditions  (air)  are:  Mm=6. 5,  TX=225.56K,  px  =  1 .44 x  1 0”1  kg  / m3 ,  rn  =  2.5mm  ,  the  cone  half 
angle  is  7  degrees,  TW=300K  (isothermal  wall),  Rew  =  19.1x1 06m~] ,  the  first  part  of  the  model 
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without  the  flare  is  1.0  meter  long.  Fig.  10  shows  our  steady  flow  solutions  computed  by  the 
fifth-order  shock  fitting  code.  The  outer  computational  boundary  is  the  bow  shock  obtained  by  the 
simulation.  The  figure  shows  a  very  good  resolution  of  the  flow  field. 


Fig.  10  Steady  flow  Mach  number  contours  for  the  FRESH  model  with  Mach  6.5  over  a  7  deg.  half  angle  blunt 
cone.  The  coordinates  are  non-dimensionalized  by  the  cone  nose  radius  of  2.5  mm. 


4  SUMMARY  OF  RESEARCH  ACCOMPLISHMENTS 

In  the  report  period,  we  have  mainly  focused  on  the  DNS  studies  of  receptivity  of  supersonic  and 
hypersonic  boundary  layer  flows.  DNS  results  are  compared  with  Maslov’s  leading-edge  receptivity 
experiment  on  Mach  5.92  flow  over  a  flat  plate,  and  Stetson's  1984  stability  experiment  on  Mach 
7.99  flow  over  a  blunt  cone.  We  also  collaborate  with  Prof.  Tumin  in  University  of  Arizona  to 
compare  numerical  and  theoretical  results  on  receptivity  of  a  Mach  8  flow  over  a  sharp  wedge  to 
wall  blowing-suction,  and  to  deeply  analyze  the  receptivity  characteristics.  Our  numerical  studies 
have  been  validated  to  be  of  high  accuracy  and  led  to  further  understanding  of  hypersonic  boundary 
layer  receptivity  mechanisms.  The  main  research  accomplishments  are 

1.  We  have  conducted  DNS  studies  on  the  receptivity  of  a  Mach  8.0  flow  over  a  sharp  wedge  of  a 
half-angle  of  5.3°  to  wall  blowing-suction.  The  results  show  that  Mode  F,  mode  S,  and  acoustic 
modes  are  simultaneously  excited  by  wall  blowing-suction  disturbances.  Downstream  of  the 
forcing  region,  mode  F  decays  owing  to  its  inherent  stability  whereas  mode  S  grows  because  of 
its  instability.  Acoustic  modes  radiate  into  the  external  flow  outside  the  boundary  layer.  As  a 
result,  mode  S  eventually  becomes  the  dominant  wave  mode  in  the  boundary  layer.  It  is  also 
found  that  the  frequency  of  the  blowing-suction  actuator  has  a  significant  effect  on  the 
receptivity  process  owing  to  the  frequency  dependence  of  the  synchronization  point  between 
mode  F  and  mode  S.  Furthermore,  the  excitation  of  mode  S  is  strongly  affected  by  profile  and 
length  of  the  forcing  actuator.  All  cases  of  numerical  simulations  consistently  show  that  the 
synchronization  point  plays  an  important  role  in  the  excitation  of  mode  S  by  wall 
blowing-suction,  i.e.,  mode  S  is  strongly  excited  only  when  the  blowing-suction  actuator  is 
located  upstream  of  the  synchronization  point.  When  the  forcing  actuator  is  downstream  of  the 
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synchronization  point,  there  is  very  little  excitation  of  mode  S,  despite  the  fact  that  the 
blowing-suction  actuator  is  still  located  within  the  unstable  region  of  mode  S. 

2.  We  collaborated  with  Prof.  Tumin  in  University  of  Arizona  to  compare  numerical  and 
theoretical  results  on  receptivity  of  the  Mach  8  flow  over  a  sharp  wedge  to  wall  blowing-suction, 
and  to  deeply  analyze  the  receptivity  characteristics.  The  perturbation  flow  field  downstream  of 
the  blowing-suction  actuator  was  decomposed  into  boundary-layer  wave  modes  with  the  help  of 
the  biorthogonal  eigenfunction  system.  It  was  shown  that  amplitudes  of  boundary-layer  wave 
modes  calculated  with  the  help  of  the  theoretical  receptivity  model  had  good  agreements  with 
those  obtained  by  projecting  the  numerical  results  onto  the  boundary-layer  wave  modes 

3.  We  have  studied  the  receptivity  of  a  Mach  5.92  over  a  flat  plate  to  various  two-dimensional  wall 
perturbations:  oscillation,  blowing-suction,  and  temperature  perturbation.  The  numerical  results 
show  that  all  types  of  two-dimensional  wall  perturbations  eventually  result  in  the  same  type  of 
instability  wave  (mode  S)  in  the  boundary  layer.  The  hypersonic  boundary-layer  flow  is  most 
sensitive  to  wall  blowing-suction  and  least  sensitive  to  wall  temperature  perturbation.  It  is  also 
noticed  that  the  profile  of  temperature  perturbation  has  a  strong  effect  on  receptivity  process. 
The  receptivity  of  the  hypersonic  boundary  layer  to  temperature  perturbation  is  proportional  to 
the  energy  introduced  to  the  steady  base  flow,  where  the  energy  is  determined  by  the  profile  of 
temperature  perturbation. 

4.  We  have  conducted  DNS  studies  on  the  transient  growth  phenomena  of  a  Mach  5.92  flat-plate 
boundary  layer  and  a  Mach  7.99  flow  over  a  blunt  cone  to  small-scale  three-dimensional 
stationaiy  roughness  element.  The  surface  roughness  is  periodic  in  spanwise  direction.  Effect  of 
thermal  flat  plate  boundary  conditions  on  the  receptivity  process  is  considered  by  comparing 
numerical  simulations  on  adiabatic  and  isothermal  flat  plates,  respectively.  The  numerical 
results  show  that  counter  rotating  streamwise  vortices  are  excited  by  stationary  roughness 
element.  The  temperature  increases  in  the  wake  region  of  roughness  peak  whereas  it  decreases 
in  the  wake  region  of  roughness  valley.  Furthermore,  the  speed  of  wake  from  roughness  peak  is 
slower  than  that  from  roughness  valley.  All  these  results  are  consistent  with  those  of  Tumin's 
theoretical  analysis.  It  is  also  concluded  based  on  the  numerical  results  that  the  roughness 
element  on  adiabatic  flat  plate  is  more  efficient  in  the  excitations  of  streamwise  vortices  and 
transient  growth. 

5.  We  developed  a  new  high-order  immersed  interface  method  for  computing  boundary  layer 
transition  with  surface  roughness.  The  work  was  motivated  by  discussions  among  the 
participants  in  the  Transition  Study  Group  Open  Forum  Held  in  the  AIAA  meeting  in  Reno  in 
January  2005.  It  was  pointed  out  that  that  there  is  a  need  to  compute  boundary  layer  stability 
and  transient  growth  with  distributive  surface  roughness.  The  new  method  can  be  arbitrarily 
high-order  accuracy  in  the  whole  flow  field,  including  the  interface  with  discontinuity,  and  it 
can  be  useful  for  the  simulations  of  transient  growth  to  finite  or  nonlinear  surface  roughness. 

5  RESULTS  OF  HYPERSONIC  BOUNDARY  LAYER  RECEPTIVITY 
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One  of  the  major  accomplishments  of  the  current  research  project  is  on  revealing  the  receptivity 
mechanism  of  a  hypersonic  boundary  layer  to  wall  blowing-suction,  and  the  effect  of 
blowing-suction  locations  of  forcing  regions.  These  results  are  discussed  in  more  detail  in  this 
section.  Details  of  these  results  have  been  published  in  Publications:  [2],  [6],  [13],  [15],  [20]. 

5.1  Flow  conditions  and  blowing-suction  model 

With  the  development  of  advanced  computers  and  numerical  techniques,  numerical  simulation  of 
the  receptivity  process  by  directly  solving  Navier-Stokes  equations  has  become  feasible.  By  solving 
the  compressible  linearized  Navier-Stokes  equations,  Malik  et  al.  [15]  investigated  the  responses  of 
a  Mach  8  flow  over  a  sharp  wedge  of  a  half-angle  of  5.3°  to  three  types  of  external  forcing:  a  planar 
free-stream  acoustic  wave,  a  narrow  acoustic  beam  enforced  on  the  bow  shock  near  the  leading 
edge,  and  a  blowing-suction  slot  on  the  wedge  surface.  They  concluded  that  these  three  types  of 
forcing  eventually  resulted  in  the  same  type  of  instability  waves  in  the  boundary  layer.  However  the 
receptivity  mechanism  was  not  studied  in  detail.  Ma  and  Zhong  [16]  studied  the  receptivity 
mechanisms  of  the  same  hypersonic  boundary  layer  to  various  freestream  disturbances,  i.e.,  fast  and 
slow  acoustic  waves,  vorticity  waves,  and  entropy  waves,  by  solving  the  two-dimensional 
compressible  Navier-Stokes  equations.  They  found  that  the  stable  modes  in  the  boundary  layer 
played  a  very  important  role  in  the  receptivity  process. 


Fig.  11  A  schematic  of  the  receptivity  of  a  hypersonic  boundary  layer  over  a  sharp  wedge  to  wall 
blowing-suction:  1.  boundary  layer;  2.  sharp  wedge;  3.  bow  shock;  4.  blowing-suction  actuator. 


In  this  project,  the  receptivity  mechanism  of  the  Mach  8  flow  over  the  sharp  wedge  to  wall 
blowing-suction  is  investigated  by  series  of  numerical  simulations.  Fig.  11  shows  a  schematic  of  the 
receptivity  of  a  hypersonic  boundary  layer  over  a  sharp  wedge  to  wall  blowing-suction.  Due  to  the 
symmetric  geometry,  only  the  flow  field  over  the  upper  surface  of  the  sharp  wedge  is  simulated. 
The  steady  base  flow  is  simulated  by  solving  the  two-dimensional  compressible  Navier-Stokes 
equations  with  a  combination  of  a  fifth-order  shock-fitting  finite  difference  method  and  a 
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second-order  TVD  scheme.  In  the  receptivity  simulations,  periodic  blowing-suction  disturbances 
are  introduced  to  the  steady  base  flow  through  a  forcing  actuator  on  the  wedge  surface.  The 
characteristics  of  boundary-layer  wave  modes  are  identified  and  evaluated  by  comparing  the  results 
of  LST  and  numerical  simulations.  The  effects  of  frequency  and  location  of  the  forcing  actuator  on 
the  receptivity  process  are  studied  by  considering  seven  cases  of  different  blowing-suction  locations 
with  1 5  frequencies  for  each  case. 

The  freestream  conditions  for  the  currently  studied  hypersonic  flow  over  the  sharp  wedge  are  the 
same  as  those  used  by  Malik  et  al.  [15],  and  Ma  and  Zhong  [16],  i.e., 


Mx  =  8.0,  Tl  =  54.78  K 

p\  =389  Pa,  Pr  =  0.72 

Reoo  ~  pl>uU mL  =8.2x106  m"' 

Total  length  of  the  wedge  surface  « 1 .5  m 


For  the  simulation  of  steady  base  flow,  the  wall  is  adiabatic,  and  the  physical  boundary  condition  of 
velocity  on  the  wedge  surface  is  the  non-slip  condition.  When  periodic  blowing-suction 
disturbances  are  enforced  on  the  steady  base  flow,  the  isothermal  temperature  condition  is  applied 
on  the  wall.  This  temperature  condition  is  a  standard  boundary  condition  for  theoretical  and 
numerical  studies  of  high  frequency  disturbances.  Meanwhile,  non-slip  condition  is  applied  on  the 
wall  except  the  forcing  region.  Inlet  conditions  are  specified.  High-order  extrapolation  is  used  for 
outlet  conditions  because  the  flow  is  supersonic  at  the  exit  boundary  except  a  small  region  near  the 
wedge  surface. 


In  Fedorov  and  Khokhlov's  [23]  theoretical  and  numerical  analysis  of  boundary-layer  receptivity  to 
wall  disturbances,  wall  blowing-suction  disturbance  has  the  following  traveling-wave  form: 
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"O' 

v' 

=  g(x) 
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w' 

0 

0'_ 
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( iacx+if}cz~ia)t ) 


(1) 


where  u' ,  v',  w'  are  velocity  disturbances  in  streamwise,  wall-normal,  and  spanwise  directions, 

respectively,  while  O'  is  the  temperature  disturbance.  The  parameters  ac  and  /?c  are  wave 

number  components  in  streamwise  and  spanwise  directions,  and  co  is  the  circular  frequency.  The 

function  g-(x)  represents  the  disturbance  profile.  Equation  (1)  implies  that  wall  blowing-suction  is 

only  related  to  the  wall-normal  velocity  disturbance.  For  this  model,  small  instantaneous  mass  flux 
is  generally  introduced  to  the  boundary  layer. 
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Another  model  of  blowing-suction  disturbances  is  to  enforce  the  mass  flux  oscillations  on  the  wall, 
which  has  been  used  by  Eibler  and  Bestek  [24],  and  Egorov  et  al.  [25].  For  this  model,  wall 
blowing-suction  is  expressed  as 


AAV 

P-K 


=  qw(x,t)  =  Asin 


X2~XlJ 


sin<y/ 


(2) 


where  A  is  the  amplitude  of  mass  flux  oscillation,  pv  and  vw  are  the  local  density  and 


wall-normal  velocity  at  the  location  x  on  the  wall,  pa  and  ux  are  the  density  and  stream-wise 


velocity  in  free-stream.  The  coordinates  x,  and  x2  represent  the  leading  and  trailing  edges  of  the 

blowing-suction  actuator.  Due  to  the  anti-symmetric  profile  of  sine  function  within  the  forcing 
region,  the  net  mass  flux  introduced  to  the  boundary  layer  is  zero  at  any  instant. 


In  current  study,  a  blowing-suction  model  similar  to  Eq.  (2)  is  used.  The  mass  flux  oscillations  on 
the  wedge  surface  within  the  blowing-suction  region  can  be  written  as: 

(p'v)'  =  qlsp(l)YJS\n(oy  (3) 

n=\ 


where  q'a  is  a  local  constant  depending  on  the  location  of  the  blowing-suction  actuator,  and  e  is  a 
small  dimensionless  parameter  representing  the  amplitude  of  the  mass  flux  oscillation.  The  function 
/?(/)  is  the  profile  function  defined  within  the  forcing  region  as 


m= 


f  20.25/5 -35.4375Z4  +15.1875/2  if  /<1 
[-20.25(2-/)’  +35.4375(2-/)4  - 1 5. 1 875(2 - /)2 


if  />  1 


The  variable  /  is  a  dimensionless  coordinate  defined  within  the  blowing-suction  region, 


/  = 


2  (s-s,) 
*  * 

Se~S, 


(4) 

(5) 


where  si ,  s'e  are  the  coordinates  of  the  leading  and  trailing  edges  of  the  blowing-suction  actuator. 

Compared  with  the  sine  profile  function  in  Eq.  (2),  the  specific  5th-order-polynomial  profile 
function  makes  smoother  the  mass  flux  oscillations  at  the  edges  of  the  forcing  actuator.  Due  to  the 
anti-symmetric  property  of  the  5th-order-polynomial  profile  function  within  the  blowing-suction 

region,  the  net  mass  flux  introduced  to  the  boundary  layer  is  zero  at  any  instant.  In  Eq.  (3),  <y*  is 


the  circular  frequency  of  multi-frequency  blowing-suction  disturbances  («  =  1, 2, — ,15) ,  which  is 
related  to  the  frequency  by, 

o)‘n  =  Inf  =  2nnf  (n  =  l,2,— ,15)  (6) 
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The  circular  frequency,  con ,  and  the  frequency,  f‘,  are  non-dimensionalized  according  to 


co.  = 


F_  = 
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°>nL 
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The  current  receptivity  studies  are  focused  on  linear  responses  of  the  boundary  layer  to  forcing 
waves.  In  the  simulations,  the  amplitude  of  mass  flux  oscillation  is  small  enough  to  preserve  the 
linear  properties  of  the  disturbances.  The  dimensionless  amplitude  coefficient,  s  in  Eq.  (3),  is 
given  as 

s  =  1.0x1 0~5  (9) 


5.2  Steady  base  flow 

The  steady  base  flow  is  simulated  by  solving  the  two-dimensional  compressible  Navier-Stokes 
equations  with  a  combination  of  a  fifth-order  shock-fitting  finite  difference  method  and  a 
second-order  TVD  scheme.  In  the  leading  edge  region,  there  exists  a  singular  point  at  the  tip  of  the 
wedge,  which  will  introduce  numerical  instability  if  the  fifth-order  shock-fitting  method  is  used  to 
simulate  the  flow.  Therefore,  the  computational  domain  for  the  fifth-order  shock-fitting  method 

starts  at  s*  =  0.00409  m  and  ends  at  s*  =  1 .48784  m ,  respectively.  In  actual  simulations,  the 

computational  domain  is  divided  into  30  zones  with  a  total  of  5936  grid  points  in  stream-wise 
direction  and  121  grid  points  in  wall-normal  direction.  Forty-one  points  are  used  in  the  overlap 
region  between  two  neighboring  zones,  which  is  proved  to  be  sufficient  to  make  the  solution 
accurate  and  smooth  within  the  whole  domain.  An  exponential  stretching  function  is  used  in  the 
wall-normal  direction  to  cluster  more  points  inside  the  boundary  layer.  On  the  other  hand,  the  grid 
points  are  uniformly  distributed  in  stream-wise  direction.  The  spatial  convergence  of  the  results 
based  on  this  grid  structure  has  been  evaluated  by  grid  refinement  studies  to  ensure  the  grid 
independence  of  the  numerical  simulations.  The  results  of  the  evaluation  of  numerical  accuracy  can 
be  found  in  Ma  and  Zhong's  paper  [16],  they  are  not  repeated  here. 

For  the  first  zone  of  the  shock-fitting  calculations,  the  inlet  conditions  are  obtained  from  the  results 
of  the  second-order  TVD  shock-capturing  scheme  which  is  used  to  simulate  the  steady  base  flow  in 
a  small  region  including  the  leading  edge.  For  other  zones,  inlet  conditions  are  interpolated  from  the 
results  of  the  previous  zone.  Fig.  12  shows  the  wall-normal  velocity  and  density  contours  near  the 
leading  edge  of  the  steady  base  flow  obtained  by  the  second-order  TVD  scheme  and  the  fifth-order 
shock-fitting  method.  The  flow  field  including  the  leading  edge  is  simulated  by  the  TVD  scheme, 

whereas  the  flow  field  after  s‘  =  0.00409  m  is  simulated  by  the  shock-fitting  method.  The  region 
between  s*  =  0.00409  m  to  s*  =  0.0064  m  is  an  overlap  region  where  the  flow  is  computed  by 
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both  methods.  This  figure  shows  that  wall-normal  velocity  and  density  contours  have  good 
agreements  within  the  overlap  region,  which  indicates  that  the  TVD  solutions  are  accurate  to  be 
used  as  inlet  conditions  for  the  fifth-order  shock-fitting  simulation  in  the  first  zone.  The 
combination  of  the  shock-fitting  method  and  the  TVD  scheme  has  also  been  validated  in 
simulations  of  supersonic  and  hypersonic  steady  base  flows  over  a  flat  plate  by  Ma  and  Zhong  [26], 
and  Wang  and  Zhong  [4], 


x*  (m)  x‘  (m) 


Fig.  12  Numerical  results  of  steady  base  flow:  1)  wall-normal  velocity  and  density  contours  near  the  leading 
edge;  2)  pressure  contour;  3)  pressure  distributions  along  the  wedge  surface;  4)  shock  front  position  and 
distribution  of  the  shock  angle;  5)  distribution  of  Mach  number  behind  the  bow  shock. 

Fig.  12  shows  the  pressure  contours  of  the  steady  base  flow  simulated  by  the  fifth-order 
shock-fitting  finite  difference  method.  The  upper  boundary  of  the  flow  field  represents  the  bow 
shock  induced  by  the  sharp  wedge  and  the  displacement  thickness  of  the  boundary  layer.  The  lower 
boundary  is  the  surface  of  the  sharp  wedge  of  a  half  angle  of  5.3°.  A  part  of  the  pressure  field  from 

x*  =  0.05  m  to  x*  =  0.15  m  is  amplified  to  show  clearly  the  pressure  contour  within  the 

boundary  layer.  It  is  noticed  that  pressure  is  approximately  a  constant  across  the  boundary  layer  and 
along  the  Mach  lines,  which  is  consistent  with  the  theories  of  the  boundary  layer  flow  and  inviscid 

supersonic  aerodynamics.  At  a  fixed  location  (constant  x*),  pressure  behind  the  shock  is  higher 

than  that  on  the  wedge  surface  due  to  the  existence  of  the  bow  shock.  Fig.  12  also  shows  the 
pressure  distribution  along  the  wedge  surface  obtained  by  the  fifth-order  shock-fitting  method.  The 
theoretical  inviscid  limit  is  also  plotted  in  the  figure  for  comparison.  Near  the  leading  edge,  there 
exists  great  pressure  gradient  for  the  numerical  result,  which  is  caused  by  the  interaction  between 
the  inviscid  outer  flow  and  the  viscous  boundary  layer.  From  upstream  to  downstream,  the 


outer-flow/boundary-layer  interaction  becomes  weaker  with  the  bow  shock  moving  away  from  the 
boundary  layer.  As  a  result,  the  pressure  approaches  the  constant  value  of  the  inviscid  limit  further 
downstream.  The  figure  also  shows  that  the  pressure  of  numerical  solution  at  the  exit  of 
computational  domain  is  greater  than  that  of  the  theoretical  inviscid  limit.  This  difference  in 
pressure  is  a  result  of  the  induced  boundary  layer  thickness  on  the  wedge  surface. 

The  shock  front  position  and  the  shock  angle  obtained  by  the  numerical  simulations  are  plotted 
together  with  those  of  the  theoretical  inviscid  limit  in  Fig.  12.  The  solid  line  represents  the  shock 
front  whereas  the  dashdot  line  is  a  straight  line  representing  the  oblique  shock  induced  by  an 
inviscid  flow  passing  the  sharp  wedge  with  the  same  free-stream  conditions.  The  figure  shows  that 
the  shock  front  of  numerical  solution  is  not  a  straight  line.  The  angle  of  the  bow  shock  decreases 
from  14.793°  near  the  leading  edge  to  1 1.307°  at  the  exit  of  the  computational  domain,  however  the 
inviscid  limit  of  the  shock  angle  is  a  constant  of  11.102°.  Fig.  12  shows  the  distribution  of  Mach 
number  behind  the  bow  shock  of  the  numerical  result.  The  theoretical  inviscid  limit  of  Mach 
number  behind  the  oblique  shock  is  also  plotted  in  the  figure  for  comparison.  It  shows  that  Mach 
number  increases  from  5.911  near  leading  edge  to  6.746  at  the  exit  of  the  computational  domain. 
For  inviscid  Mach  8  flow  over  the  sharp  wedge  of  a  half-angle  of  5.3°,  Mach  number  behind  the 
shock  is  6.798.  Figure  6  shows  that  the  shock  angle  and  Mach  number  behind  the  shock  approach 
corresponding  values  of  the  theoretical  inviscid  limit  further  downstream,  where  the  interaction 
between  the  inviscid  outer  flow  and  the  viscous  boundary  layer  becomes  weak.  Again,  large 
gradients  of  the  shock  angle  and  Mach  number  behind  the  shock  near  the  leading  edge  are  caused 
by  the  interaction  between  the  inviscid  outer  flow  and  the  viscous  boundary  layer. 


5.3  Characteristics  of  boundary-layer  wave  inodes 


The  characteristics  of  boundary-layer  wave  modes  of  the  Mach  8  flow  over  the  sharp  wedge  of  a 
half-angle  of  5.3°  is  studied  by  LST  based  on  a  multi-domain  spectral  method  of  Malik  [27].  The 
LST  code  was  written  by  Ma  and  Zhong  [28].  The  velocity,  pressure,  and  temperature  disturbances 
are  represented  by  harmonic  waves  of  the  form 
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Similar  to  Eq.  (1),  the  parameters  ac  and  pc  are  wave  number  components  in  streamwise  and 


spanwise  directions,  and  coc  is  the  circular  frequency.  For  the  current  two-dimensional  flow  cases. 


Pc  =  0 .  Substituting  disturbances  of  Eq.  (10)  and  the  steady  base  flow  into  the  compressible 
linearized  Navier-Stokes  equations,  an  ordinary-differential-equation  system  is  achieved,  i.e., 
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(11) 


d2  d  \ 

A—r  +  B  —  +  C  U  =  0 
dy 2  dy  Y 

where  <j>  is  the  disturbance  vector  defined  by  jw,v,p, T’.wj  .  The  coefficient  matrices  A,  B,  and 

C  were  given  in  Malik's  paper  [27],  In  spatial  stability  analysis,  the  two  parameters,  coc  and  J3C, 

are  specified  as  real  numbers.  The  streamwise  wave  number,  ac ,  is  a  complex  number  and  solved 

as  the  eigenvalue  of  the  ordinary-differential-equation  system.  The  complex  wave  number  ac  can 
be  expressed  as 


ac=ar+  iai  (12) 

where  ai  is  the  local  growth  rate.  A  boundary-layer  wave  mode  is  unstable  when  a.  <  0  whereas 

it  is  stable  when  ai  >  0 .  The  real  part,  ar ,  is  the  local  wave  number  which  can  be  used  to  define 
the  local  phase  velocity: 
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(13) 


Both  the  wave  number  and  the  phase  velocity  can  be  used  to  identify  the  boundary-layer  wave 
mode. 


Fig.  13  Characteristics  of  boundary-layer  modes:  1)  distribution  of  dimensionless  phase  velocity:  2) 
eigenfunctions  of  boundary-layer  mode  at  synchronization  point;  3)  distribution  of  growth  rate. 

Fig.  13  shows  the  dimensionless  phase  velocities  of  boundary-layer  wave  modes  at  three 
dimensionless  frequencies  as  a  function  of  the  dimensionless  circular  frequency  co .  These  three 

frequencies  are  the  same  as  Fs ,  F6 ,  and  Fn  in  Table  2.  The  three  dashed  lines  represent  the 
dimensionless  phase  velocities  of  fast  acoustic  mode  (a  =  l  +  M^'),  vorticity  or  entropy  mode 
(a  =  l),  and  slow  acoustic  mode  (a  =  l-M^‘),  respectively.  The  excellent  agreement  of  the  phase 
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velocities  at  different  dimensionless  frequencies  indicates  that  the  phase  velocity  is  a  function  of  co 
only.  The  figure  clearly  shows  that  mode  F  originates  from  fast  acoustic  mode.  As  co  increases, 
the  phase  velocity  of  mode  F  decreases  and  intersects  that  of  vorticity  or  entropy  mode.  On  the 
other  hand,  mode  S  originates  from  slow  acoustic  mode.  The  figure  also  shows  that  mode  S 

synchronizes  with  mode  F  at  the  point  of  6^=0.11443  and  as  =  0.93349 .  At  the  synchronization 

point,  the  dimensionless  phase  velocities  of  mode  S  and  mode  F  are  the  same,  and  their 
eigenfunctions  have  similar  profile.  The  figure  shows  that  the  flow  fluctuations  are  confined  within 
the  boundary  layer.  Except  the  differences  near  the  edge  of  the  boundary  layer,  the  eigenfunctions 
of  mode  F  and  mode  S  agree  very  well  at  the  synchronization  point.  Although  Fig.  13  shows  that 

the  synchronization  point  at  different  frequencies  have  a  constant  value  of  6^=  0.1 1443) ,  the 
dimensional  location  of  the  synchronization  point,  s]n,  are  different  for  different  dimensionless 


frequencies.  The  synchronization  location  in  s'  coordinate  for  a  given  dimensionless  frequency 


can  be  calculated  using  the  following  formula 


Mf.y 


Fig.  13  also  shows  the  growth  rates  of  boundary-layer  wave  modes  at  the  same  set  of  three 
frequencies  as  a  function  of  the  dimensionless  circular  frequency.  The  vertical  dashed  line 
represents  the  position  of  the  synchronization  point,  whereas  the  horizontal  dashdot  line  stands  for 

the  neutral  modes  («,  =0).  The  growth  rates  of  either  mode  S  or  mode  F  of  the  three  frequencies 

are  almost  the  same,  and  they  are  approximately  functions  of  co  only.  Mode  S  is  unstable  in  the 

region  from  co,  =  0.04  to  co„  =  0.23  whereas  mode  F  is  always  stable.  For  unstable  mode  S,  it  is 

also  noticed  that  the  growth  rate  upstream  of  the  synchronization  point  (co<  0.095)  is  much  larger 
than  those  around  and  downstream  of  the  synchronization  point.  This  indicates  that  mode  S  is  more 

unstable  around  and  downstream  of  the  synchronization  point.  The  parameters  co,  and  co„  are 

called  the  Branch  I  and  Branch  II  neutral  points  of  mode  S.  Mode  S  is  stable  upstream  of  Branch  I 
neutral  point  and  downstream  of  Branch  II  neutral  point.  The  location  of  the  Branch  II  neutral  point 

in  s'  coordinate,  which  changes  with  different  dimensionless  frequencies,  can  be  calculated  by 


Equations  (14)  and  (15)  show  that  when  F  increases,  the  corresponding  s]n  and  s'„n  decrease.  In 

other  words,  the  synchronization  point  and  the  Branch  II  neutral  point  move  upstream  when  F 
increases.  Table  2  lists  dimensional  frequency,  dimensionless  circular  frequency,  dimensionless 
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frequency,  locations  of  synchronization  point  and  Branch  II  neutral  point  for  the  15  sets  of 
frequencies  considered  in  the  current  study 


Table  2  Dimensional  frequency,  dimensionless  circular  frequency.,  dimensionless  frequency, 
locations  of  synchronization  point  and  Branch  II  neutral  point  for  the  15  sets  of  frequencies 

considered  in  the  current  study 


H 

/>  Hz) 

<4  (kHz) 

F„xl06 

4„(m) 

4„(m) 

n 

14.92 

93.74 

9.63 

16.9885 

68.38297 

2 

29.84 

187.48 

19.26 

4.2471 

17.09574 

3 

44.76 

28.89 

1.8876 

7.59811 

4 

59.68 

374.97 

38.52 

1.0618 

4.27394 

5 

74.60 

48.15 

0.6795 

2.73532 

6 

89.52 

| 

57.78 

0.4719 

1.89953 

ma 

104.44 

656.19 

67.41 

0.3467 

1.39557 

8 

119.36 

749.94 

77.04 

0.2654 

1.06848 

9 

134.28 

843.68 

86.67 

0.2097 

0.84423 

B9 

149.20 

937.42 

96.30 

0.1699 

0.68383 

TO 

164.12 

1031.16 

105.93 

0.1404 

0.56515 

TO 

179.04 

1124.91 

115.56 

0.1180 

0.47488 

13 

193.96 

1218.65 

125.19 

0.1005 

0.40463 

14 

1312.39 

134.82 

0.0867 

0.34889 

15 

223.80 

1406.13 

144.45 

0.0755 

0.30392 

5.4  Receptivity  to  a  single-frequency  blowing-suction  actuator 

The  receptivity  result  of  the  hypersonic  boundary  layer  to  a  blowing-suction  actuator  at  a  single 
frequency  is  firstly  presented  in  this  section.  The  forcing  disturbance  at  the  frequency  of 

fl  =  74.60  kHz  is  introduced  through  a  blowing-suction  actuator  on  the  wedge  surface  from 
s’  =0.10184  m  to  s’  =0.11384  m.  Table  3  lists  the  parameters  of  the  blowing-suction  actuators 

at  seven  different  locations.  The  dimensional  values  of  q’0  and  other  parameters  of  the  forcing 

actuator  considered  in  this  section  are  the  same  as  case  2  in  Table  3.  The  value  of  e  is  given  in  Eq. 
(9).  Due  to  the  fact  that  only  very  weak  disturbances  are  considered,  the  current  receptivity  is  in  the 
linear  region. 

Fig.  14  shows  the  distribution  of  dimensionless  instantaneous  pressure  perturbation  along  the  wedge 
surface.  The  amplification  of  the  pressure  perturbation  from  upstream  to  downstream  indicates  the 
excitation  of  unstable  modes  in  the  boundary  layer.  The  significant  growth  of  the  pressure 
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perturbation  starts  at  a  location  close  to  the  corresponding  synchronization  point  at  s'sS  =  0.6795  m. 

It  is  very  likely  that  mode  S  is  excited  in  the  boundary  layer  by  the  blowing-suction  disturbance, 
because  mode  S  is  more  unstable  around  and  downstream  of  the  synchronization  point. 

Table  3  Constant  Error!  Objects  cannot  be  created  from  editing  field  codes,  and  other  parameters 
of  blowing-suction  actuator  for  the  seven  cases  by  which  the  effect  of  location  of  the 
blowing-suction  actuator  on  receptivity  is  investigated 


case 

%  (kg/m2s) 

se  (m) 

s’(m> 

1 

0.214139 

0.05184 

0.06384 

0.05784 

2 

0.125188 

0.10184 

0.11384 

0.10784 

3 

0.096130 

0.15184 

0.16384 

0.15784 

4 

0.080666 

0.20184 

0.21384 

0.20784 

5 

0.070759 

0.25184 

0.26384 

0.25784 

6 

0.063745 

0.30184 

0.31384 

0.30784 

7 

0.054453 

0.40184 

0.41384 

0.40784 

tL 

II 

0.0002 

:  I 

00001 
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Fig.  14  Distribution  of  dimensionless  pressure  perturbation  along  the  wedge  surface  and  comparisons  of 
numerical  and  LST  local  wave  number  and  local  growth  rate. 

In  order  to  check  the  properties  of  the  unstable  modes,  a  Fourier  transform  is  applied  to  the 
dimensionless  instantaneous  pressure  perturbation  along  the  wedge  surface,  which  leads  to 

pX*\t)  =  Z|p,.(S-)K[',-,'V'y]  (16) 

where  p'(s\t *)  represents  the  dimensionless  instantaneous  pressure  perturbation  along  the  wedge 


surface.  In  above  equation,  |/?'n(s‘)|  and  0’n(s‘)  are  the  perturbation  amplitude  and  phase  angle, 
respectively.  Once  |/?'n(/)|  and  are  known,  a  local  wave  number  ( am )  and  a  local 


growth  rate  ( am )  of  the  perturbation  at  the  frequency  /„*  can  be  defined  as 
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(17) 


L'  d\p'„  1 
|  P'n\  ds' 


(18) 


where  Z*  is  the  length  scale  of  local  boundary  layer  thickness  as  defined  by  Eq.  (7).  The 
parameters  am  and  am  will  represent  the  true  wave  number  and  growth  rate  only  if  the 


disturbance  is  dominated  by  a  single  wave  mode.  Otherwise,  the  disurbance  needs  to  be 
decomposed  in  order  to  check  properties  of  a  specific  mode.  For  example,  Tumin,  Wang,  and 
Zhong  [2]  decomposed  the  perturbation  at  a  location  just  downstream  of  the  blowing-suction 
actuator  with  a  biorthogonal  eigenfunction  system,  where  mode  F,  mode  S,  and  acoustic  modes 
coexisted  and  none  of  them  was  dominant. 


Since  our  focus  is  on  numerical  simulation  results,  we  only  consider  the  later  stage  where  mode  S  is 
the  dominant  mode  in  the  boundary  layer.  In  this  case,  Eqs.  (17)  and  (18)  can  be  used  to  check  the 
properties  of  the  unstable  mode  S.  Fig.  14  compares  local  wave  numbers  and  local  growth  rates  at 

the  frequency  /5’  obtained  by  Eqs.  (17)  and  (18)  with  the  corresponding  values  of  mode  F  and 

mode  S  computed  by  LST,  respectively.  These  two  figures  show  that  wave  numbers  and  growth 
rates  of  the  numerical  results  agree  well  with  those  of  mode  S  after  the  location  of  o  =  0.095 

(/  =  0.477  m).  As  shown  in  Fig.  13,  this  is  the  location  where  the  growth  rate  of  mode  S  decrease 

significantly  and  mode  S  grows  to  the  dominant  mode  in  the  boundary  layer,  which  indicates  that 
the  unstable  mode  excited  by  the  blowing-suction  forcing  is  mode  S.  It  is  also  shown  that  the  two 
sets  of  growth  rates  agree  very  well  for  co  from  0.102  to  0.131.  When  co  is  larger  than  0.131,  the 
growth  rate  from  numerical  simulation  is  larger  than  that  from  LST,  which  means  that  mode  S 
obtained  by  numerical  simulation  becomes  more  stable  than  that  predicted  by  LST.  In  other  words, 
the  Branch  II  neutral  point  of  mode  S  obtained  by  numerical  simulation  moves  upstream  compared 
with  that  predicted  by  LST.  A  possible  reason  is  due  to  the  non-parallel  effects  that  are  not 
considered  in  linear  stability  theory.  This  result  is  consistent  for  all  cases  considered  in  the  current 
study. 

Fig.  15  compares  the  eigenfunctions  of  mode  F  and  mode  S  with  the  amplitude  of  pressure 
perturbation  obtained  from  the  numerical  simulation  at  four  different  locations.  Figure  (a)  shows 
that  the  numerical  pressure  perturbation  is  not  confined  within  the  boundary  layer,  and  it  is  quite 
different  from  the  eigenfunctions  of  mode  F  and  mode  S.  These  characteristics  of  pressure 
perturbation  are  caused  by  the  coexistence  of  mode  F,  mode  S,  and  acoustic  modes  in  the  boundary 
layer  just  downstream  of  the  blowing-suction  region.  Figure  (b)  shows  that  the  numerical  pressure 
perturbation  is  much  closer  to  the  eigenfunction  of  mode  S,  and  it  is  almost  confined  within  the 
boundary  layer.  These  characteristics  of  pressure  perturbation  indicate  that  mode  S  becomes  the 
dominant  mode  in  the  local  boundary  layer  with  acoustic  modes  radiating  into  the  external  flow 
outside  the  boundary  layer  and  mode  F  decaying  owing  to  its  inherent  stability.  Figures  (c)  and  (d) 
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show  good  agreements  between  the  eigenfunctions  of  mode  S  and  the  pressure  perturbation 
obtained  from  the  numerical  simulation,  because  mode  S  is  the  dominant  mode  in  the  local 

boundary  layer  downstream  of  s’  ~  0.6  m. 


Ip'I'IPo'I  Ip'KIPo'I 

(c)  (d) 

Fig.  15  Comparisons  of  the  eigenfunctions  of  mode  F  and  mode  S  with  the  amplitude  of  pressure  perturbation 
obtained  from  numerical  simulation  at  different  locations:  (a)  ry=0.070872  (s*=0.26534  m);  (b) 

<y=0.093855  (s’ =0.46534  m);  (c)  co= 0.11223  (s*  =0.66534  m);  (d)  eo= 0.14530  (s*=l. 11534  m). 

The  numerical  results  and  analysis  above  show  that  mode  F,  mode  S,  and  acoustic  modes  are 
simultaneously  excited  by  the  blowing-suction  disturbances.  All  these  modes  coexist  in  the 
boundary  layer  just  downstream  of  the  forcing  region,  which  leads  to  the  strong  modulations  of 
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pressure  perturbation  amplitude.  Far  downstream  of  the  forcing  region,  acoustic  modes  radiate  into 
the  external  flow  outside  the  boundary  layer.  Mode  F  decays  owing  to  its  inherent  stability  whereas 
mode  S  grows  substantially  because  of  its  instability.  As  a  result,  mode  S  becomes  the  dominant 
mode  in  the  boundary  layer. 

5.5  Effect  of  location  of  the  blowing-suction  actuator 

In  order  to  investigate  the  effect  of  location  of  the  blowing-suction  actuator  on  receptivity,  a  series 
of  numerical  simulations  have  been  carried  out  for  different  actuator  locations.  Specifically,  seven 
cases  of  different  actuator  locations  are  considered.  Fig.  16  shows  a  schematic  of  forcing  actuator 

locations  for  the  seven  cases.  In  each  case,  wall  blowing-suction  with  fifteen  frequencies  ( f*  listed 

in  Table  2  with  n  from  1  to  15)  is  introduced  on  the  wedge  surface.  The  subsequent  responses  of  the 
boundary  layer  are  simulated  by  the  fifth-order  shock-fitting  finite  difference  method. 
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Fig.  16  A  schematic  of  locations  of  the  blowing-suction  actuator  for  the  seven  cases:  1 .  bow  shock;  2. 
boundary  layer;  3.  blowing-suction  actuator;  4.  sharp  wedge. 


The  model  of  the  multi-frequency  blowing-suction  actuator  is  given  by  Eq.  (3).  The  constant  q'0  in 


model  equation  are  different  for  the  seven  cases  because  it  is  locally  defined  at  the  leading  edge  of 
the  blowing-suction  actuator.  The  location  of  the  forcing  actuator  is  defined  as 


s.  +  S„ 


s.  =  - 


(19) 


Constant  q’0  and  other  parameters  of  the  blowing-suction  actuator  for  the  seven  cases  considered  in 
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the  current  section  are  listed  in  Table  3.  As  already  been  mentioned,  $qA*_0$  is  defined  at  the 
location  s’  as  the  multiple  of  unperturbed  density  on  the  wall  and  wall-normal  velocity  after  the 
bow  shock  with  the  dimension  of  kg/m2s. 

It  is  noticed  that  the  relative  location  of  the  blowing-suction  actuator  with  respect  to  the 
synchronization  point  plays  a  very  important  role  in  the  receptivity  process.  According  to  Fig.  13, 
the  synchronization  point  between  mode  F  and  mode  S  has  a  dimensionless  circular  frequency  of 

<^=0.11443.  Although  the  synchronization  point  in  a>  coordinate  does  not  depend  on  the 
dimensional  frequency,  its  location  in  s’  coordinate  is  different  for  different  dimensional 
frequency.  The  synchronization  points  in  s’  coordinate  for  the  1 5  frequencies  are  calculated  by  Eq. 
(14)  and  tabulated  in  Table  2  as  s’m . 


Table  4  Relative  location  of  the  blowing-suction  actuator  with  respect  to  the  synchronization 
point  for  the  seven  cases  by  which  the  effect  of  location  of  the  blowing-suction  actuator  on 

receptivity  is  investigated 


case 

•  * 

s  <s 

c  srt 

•  * 

s  >  s 

c  srt 

1 

n  from  1  to  15 

none 

2 

n  from  1  to  12 

n  from  13  to  15 

3 

n  from  1  to  10 

n  from  11  to  15 

4 

n  from  1  to  9 

n  from  10  to  15 

5 

n  from  1  to  8 

n  from  9  to  15 

6 

n  from  1  to  7 

n  from  8  to  15 

7 

n  from  1  to  6 

n  from  7  to  15 

Table  2  shows  that  the  synchronization  point  moves  upstream  when  the  dimensional  frequency 
increases.  Therefore,  the  blowing-suction  actuators  for  the  seven  cases  are  located  either  upstream 

or  downstream  of  the  corresponding  synchronization  point  of  f* .  Table  4  lists  relative  locations  of 

the  blowing-suction  actuators  with  respect  to  the  synchronization  points  for  the  seven  cases. 

“s’  <  s’m  ”  represents  that  the  forcing  actuator  is  located  upstream  of  the  synchronization  point  of 

f*  ,  while  “  s’c  >  s’sn  ”  represents  that  the  biowing-suction  actuator  is  downstream  of  the 

synchronization  point  of  f’ .  For  example,  for  the  dimensional  frequency  of  f’  =  74.60  kHz, 
locations  of  the  blowing-suction  actuators  in  all  cases  are  upstream  of  the  synchronization  point. 


30 


While  for  the  frequency  of  =149.20  kHz,  locations  of  the  blowing-suction  actuators  in  cases 

from  1  to  3  are  upstream  of  the  synchronization  point,  whereas  locations  of  the  blowing-suction 
actuators  in  other  cases  are  downstream  of  the  synchronization  point. 

To  show  the  effect  of  location  of  the  blowing-suction  actuator  on  receptivity  more  efficiently,  the 
perturbations  of  the  same  frequency  are  plotted  together  for  the  seven  cases.  Fig.  17  shows  the 
pressure  perturbation  amplitudes  of  the  same  frequency  for  the  seven  cases  of  different 
blowing-suction  locations.  In  these  figures,  the  numbers  of  1  to  7  within  the  small  rectangular 

represent  the  case  number.  Pressure  perturbation  amplitudes  at  the  frequencies  f‘,  /2*,  and  /3* 

are  very  small,  because  mode  S  is  slightly  unstable  in  the  region  upstream  of  the  synchronization 
point.  Therefore,  pressure  perturbations  of  these  three  frequencies  are  left  out.  For  the  frequency  of 

/4*=  59.68  kHz,  the  blowing-suction  actuators  for  all  seven  cases  are  located  upstream  of  the 

synchronization  point  at  s’4  =1.0618  m.  Figure  (a)  shows  that  mode  S  is  strongly  excited  for  all 

seven  cases.  Furthermore,  the  amplitudes  of  pressure  perturbations  decrease  significantly  when  the 
actuator  shifts  from  upstream  to  downstream.  Figures  (b)  and  (c)  show  similar  results  for  pressure 

perturbations  at  the  frequencies  of  f*=  74.60  kHz  and  /6*=  89.52  kHz,  respectively.  At  these 

two  frequencies,  Table  4  shows  that  the  blowing-suction  locations  are  all  upstream  of  the 
corresponding  synchronization  points.  Therefore,  mode  S  is  strongly  excited  for  all  seven  cases. 


When  the  frequency  is  fj=  104.44  kHz,  the  location  of  the  blowing-suction  actuator  is  upstream 

of  the  synchronization  point  at  .s’7  =  0.3467  m  for  cases  from  1  to  6,  however  it  is  downstream  of 

the  synchronization  point  for  case  7.  Figure  (d)  shows  that  mode  S  is  strongly  excited  for  cases  from 
1  to  6.  In  case  7,  there  is  very  little  excitation  of  mode  S,  despite  the  fact  that  the  blowing-suction 
actuator  is  still  located  within  the  unstable  region  of  mode  S.  The  results  of  all  cases  tested  in  this 
study  show  a  consistent  trend  that  mode  S  is  very  weakly  excited  when  blowing-suction  location  is 

downstream  of  the  synchronization  point.  For  the  frequency  of  /g*=  1 19.36  kHz,  figure  (e)  shows 

that  mode  S  is  strongly  excited  for  cases  from  1  to  5,  while  there  is  not  much  excitation  of  mode  S 
for  case  6  and  case  7.  Again,  this  is  resulted  from  the  fact  that  the  locations  of  the  blowing-suction 

actuators  for  cases  from  1  to  5  are  upstream  of  the  synchronization  point  at  s*8  =  0.2654  m, 

whereas  the  forcing  actuator  is  located  downstream  of  the  synchronization  point  for  case  6  and  case 
7.  It  needs  to  be  emphasized  that  the  blowing-suction  actuator  is  still  located  within  the  unstable 
region  of  mode  S  even  in  cases  6  and  7. 
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As  the  frequency  changes  to  /9*=  134.28  kHz  in  figure  (f),  the  blowing-suction  actuators  are 

located  upstream  of  the  synchronization  point  at  s]9  =  0.2097  m  for  cases  from  1  to  4,  whereas 

locations  of  the  forcing  actuator  for  cases  from  5  to  7  are  downstream  of  the  synchronization  point. 
The  results  for  this  frequency  are  again  consistent  with  those  of  the  previous  frequencies  regarding 
to  the  excitation  of  mode  S.  Figure  (f)  shows  that  mode  S  is  strongly  excited  for  cases  from  1  to  4, 

however  the  excited  mode  S  is  very  weak  for  cases  from  5  to  7.  For  the  frequencies  f*  with  n 

from  10  to  11,  similar  conclusion  can  be  drawn.  For  still  higher  frequencies,  the  last  three  figures 
show  that  pressure  perturbation  is  amplified  in  a  very  small  domain  downstream  of  the 
blowing-suction  actuator,  because  the  unstable  region  of  mode  S  is  small  for  higher  frequencies. 
Consistently,  mode  S  is  significantly  excited  when  the  blowing-suction  actuator  is  located  upstream 
of  the  synchronization  point.  When  the  blowing-suction  actuator  is  downstream  of  the 
synchronization  point,  there  is  a  very  little  excitation  of  mode  S. 

To  summarize,  the  current  simulation  results  have  consistently  indicated  that  the  synchronization 
point  between  mode  F  and  mode  S  plays  an  important  role  in  the  excitation  of  mode  S  by  the 
blowing-suction  actuator.  Mode  S  is  strongly  excited  only  when  the  blowing-suction  actuator  is 
located  upstream  of  the  synchronization  point.  On  the  other  hand,  when  the  forcing  actuator  is 
downstream  of  the  synchronization  point,  there  is  a  very  weak  excitation  of  mode  S.  This  happens 
even  when  the  blowing-suction  actuator  is  still  within  the  unstable  region  of  mode  S.  Therefore,  the 
synchronization  point  is  critical  to  the  receptivity  process.  The  relationship  between  the  location  of 
the  blowing-suction  actuator  and  the  synchronization  point  suggests  that,  in  order  to  control  or  delay 
the  laminar-turbulent  transition  more  efficiently,  the  blowing-suction  actuator  should  be  placed 
upstream  of  the  synchronization  point  between  mode  S  and  mode  F. 

6  DEVELOPMENT  OF  A  NEW  HIGH-ORDER  IMMERSIED  INTERFACE 
METHOD 


Another  major  accomplishment  of  the  current  research  project  is  the  development  of  a  new 
high-order  immersed  interface  method.  The  work  was  motivated  by  discussions  among  the 
participants  in  the  Transition  Study  Group  Open  Forum  Held  in  the  AIAA  meeting  in  Reno  in 
January  2005.  It  was  pointed  out  that  that  there  is  a  need  to  compute  boundary  layer  stability  and 
transient  growth  with  distributive  surface  roughness.  The  new  method  can  be  arbitrarily  high-order 
accuracy  in  the  whole  flow  field,  including  the  interface  with  discontinuity,  and  it  can  be  useful  for 
the  simulations  of  transient  growth  to  finite  or  nonlinear  surface  roughness.  The  new  high-order 
immersed  interface  method  has  been  tested  and  published  in  Publications:  [2],  [12]. 

6.1  Introduction 
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Recently,  there  has  been  strong  interest  in  developing  numerical  methods  for  computing  multi-phase 
flow  with  unsteady  interface.  These  methods  have  many  practical  applications,  such  as  the 
simulation  of  the  dynamics  of  gas  bubbles  in  a  liquid  [29],  drop  deformation  and  breakup  in  viscous 
flow  [30],  free  surface  flow  [31,  32],  and  the  breakup  of  a  liquid  jet  emanating  into  another  fluid 
[33]. 

Compared  with  single-phase  numerical  methods,  algorithms  for  two-phase  flow  simulation  face 
additional  difficulties  related  to  the  interface  treatment.  Firstly,  the  shape  of  the  interface  can  be 
complex,  and  can  undergo  change,  merge  and  breakup  during  the  course  of  the  simulation. 
Consequently,  it  is  difficult  to  use  body-fitted  unsteady  grid  to  fit  the  evolving  interface.  A  fixed 
Cartesian  grid,  where  the  interface  can  cut  through  the  grid  lines,  is  often  used.  In  a  fixed  grid,  the 
interface  can  be  treated  by,  among  others,  the  volume-of-fluid  method,  the  front  tracking  method 
[34,  35],  the  level-set  method  [36-39],  and  the  boundary  element  method  [40],  Secondly,  flow 
variables  and  their  derivatives  can  be  discontinuous  across  the  interface.  Specific  jump  conditions  at 
the  interface  depend  on  the  physical  property  of  the  problems,  the  unsteadiness  of  the  interface,  and 
the  geometric  characteristics  of  the  interface.  Consequently,  special  treatments  are  necessary  for 
computing  flow  equations  at  grid  points  adjacent  to  the  interface  (i.e.  irregular  points).  One  of  the 
popular  methods  in  treating  interface  discontinuity  is  the  immersed  boundary  method  (IBM) 
originally  developed  by  Peskin  [21]  for  simulating  blood  flow  in  the  heart.  The  basic  idea  of  the 
immersed  boundary  method  is  to  model  the  interface  by  adding  a  delta-function  source  term  to  the 
Navier-Stokes  equations.  The  resulting  equations  are  then  discretized  by  a  standard  finite  difference 
method  in  a  fixed  Cartesian  (or  non-Cartesian)  grid.  The  singular  delta  function  is  regularized  by  an 
approximate  smooth  function  spanning  a  few  grid  cells. 

An  alternative  to  the  immersed  boundary  method  is  the  “sharp-interface”  methods  which  achieve 
uniformly  second  (or  higher)  order  accuracy  by  incorporating  the  jump  conditions  into  the  finite 
difference  formulas.  The  immersed  interface  method  (IIM)  introduced  by  LeVeque  and  Li  [22]  is 
one  of  these  methods.  In  presenting  their  original  IIM  method,  LeVeque  and  Li  (1994)  considered 
finite  difference  methods  for  the  following  elliptic  equation: 

V  •  (J3(x  )Vk(x))  +  k(x)u(x)  =  f(x)  (20) 

The  equation  is  defined  in  a  simple  region  with  a  uniform  Cartesian  grid.  Fig.  1 8  shows  a  schematic 
of  a  two-dimensional  grid.  There  is  an  irregular  surface  T ,  which  may  cut  across  the  grid  lines,  in 

the  computational  domain.  Across  the  interface,  /?,  k,  and  /  may  be  discontinuous,  and  along  it 

/  may  have  a  delta  function  singularity.  In  the  derivation  of  finite  difference  formulas,  the 

computational  grid  points  are  classified  into  two  categories  depending  on  their  relative  locations 
with  respect  to  T :  regular  points  away  from  T  and  irregular  points  adjacent  to  T .  A  globally 

0(h2)  accuracy  is  achieved  by  using  the  conventional  0(h2)  central  scheme  for  the  regular 
points  and  a  locally  0(h)  scheme  for  the  irregular  points.  In  the  one-dimensional  case,  a  finite 
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difference  formula  of  0(h)  accuracy  at  an  irregular  point  uses  a  three-point  grid  stencil  together 

with  an  additional  correction  term.  A  Taylor  series  expansion  at  the  interface  is  used  to  obtain  a  set 
of  linear  equations  for  the  undetermined  coefficients  and  the  correction  term.  The  linear  equations 
are  often  problem  dependent,  and  they  need  to  be  solved  numerically  every  time  they  are  used  in  the 

simulation.  In  order  to  reach  a  locally  0(h)  approximation,  the  correction  term  requires  jump 
conditions  of  up  to  the  second  derivatives,  i.e. 

[«],  [/?wj,  and  [/?«„]  (21) 


where  [  ]  denotes  the  jump  in  variables  across  the  interface. 


Fig.  18  A  schematic  of  two-dimensional  uniform  grid  with  an  immersed  interface  T  with  discontinuity  in 

solutions. 

The  immersed  interface  methods  have  been  applied  to  the  Stokes  flow  with  elastic  boundaries  or 
surface  tension,  Hele-Shaw  flow,  incompressible  flow  based  on  the  Navier-Stokes  equations  with 
singular  source  terms,  and  nonlinear  problems  in  magneto-rheological  fluids.  Despite  these 
applications,  the  immersed  interface  methods  are  often  difficult  to  apply  to  complex  two  or  three 
dimensional  two-phase  flow  problems.  In  order  to  maintain  a  second-order  accuracy,  it  is  necessary 
to  obtain  jump  conditions  at  the  interface  for  flow  variables  and  their  first  and  second  derivatives. 
For  the  Navier-Stokes  equations  with  an  interface  of  discontinuity,  it  is  easy  to  derive  the  physical 
jump  conditions  for  flow  variables  and  their  first  derivatives  across  the  interface.  But  it  is  difficult 
to  obtain  jump  conditions  for  the  second  or  high  order  derivatives.  In  order  to  develop  third  or 
higher  order  immersed  interface  methods,  it  is  necessary  to  obtain  jump  conditions  for  the  third  and 
higher  derivatives.  In  addition,  the  finite  difference  formulas  of  the  original  immersed  interface 
method  need  to  be  re-derived  for  different  problems.  The  coefficients  and  the  correction  terms  in  the 
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finite  difference  formulas  at  irregular  points  cannot  be  obtained  explicitly.  They  are  often  computed 
numerically  by  solving  a  matrix  equation.  The  repeated  computations  for  the  coefficients  and 
correction  terms  can  be  computationally  expensive. 

Based  on  the  brief  review  above,  it  is  desirable  for  a  high-order  immersed  interface  method  to  have 
the  following  properties: 

1.  Only  two  physical  jump  conditions  of  the  variables  and  their  first  derivatives  should  be 
needed  in  the  second  or  higher  order  immersed  interface  methods. 

2.  Finite  difference  formulas  at  irregular  points  should  be  expressed  in  a  general  explicit  form 
(without  the  need  to  compute  matrix  equations  repeatedly)  so  that  they  can  be  applied  to 
different  problems  without  any  modification. 

The  derivation,  analysis,  and  test  results  of  the  new  methods  are  presented  in  following  sections,  the 
new  high-order  immersed  interface  methods  are  presented  in  this  paper  for  elliptic  equations  in  the 
form  of  Eq.  (20)  with  imbedded  interface  of  discontinuity  only.  Nevertheless,  the  method  has 
potential  advantages  in  the  application  to  two-phase  flow  because  of  its  high-order  accuracy  and 
simplicity  in  applications  by  requiring  only  the  physical  jump  conditions  for  variables  and  their  first 
derivatives  are  needed  in  the  finite  difference  formulas.  The  derivation  of  jump  conditions  involving 
the  second  or  higher  order  derivatives  can  be  difficult  for  two-phase  flow  problem  involving  the 
Navier-Stokes  equations. 

6.2  Explicit  Finite  Difference  Formulas  at  Irregular  Grid  Points 


c  h  Interface  T 


i-n+1  •••  i-1  i  i+1  •••  i+m 


h 


(  n  points  on  left  side) 


(  m  points  on  right  side) 


Fig.  19  Unifrom  grid  with  an  interface  located  on  the  right  side  of  irregular  point  i  with  a  general  grid  stencil 
of  n  and  m  points  on  the  left  and  right  sides,  respectively  . 

The  new  high-order  immersed  interface  method  is  presented  for  one-dimensional  differential  elliptic 
equations  in  the  form  of  Eq.  (20)  in  this  section.  The  method  is  extended  to  two-dimensional  elliptic 

equations  afterward.  For  simplicity,  only  the  finite  difference  approximation  to  (du / dx):  and 


36 


(d2u/dx2)i  is  presented,  though  formulas  for  higher  derivatives  can  be  easily  obtained  by  the  same 

method.  A  uniform  grid  of  mesh  size  h  shown  in  Fig.  19  is  used  for  the  discretization.  Without 
losing  generality,  it  is  assumed  that  the  origin  of  the  coordinate  system  is  located  at  grid  point 
i.e. 

xl+k=kh  (k=0,  ±1,  ±2,  •••)  (22) 


The  interface  is  located  at: 


xr  =  x.  +crh  =  crh 

where  a  is  the  interface  location  parameter  which  satisfies 

0  <  cr  <  1 


(23) 

(24) 


As  discussed  in  the  preceding  section,  only  two  jump  conditions  involving  u  and  ux  are  used  in 

the  finite  difference  approximation  of  the  derivatives.  A  general  jump  conditions  across  the  interface 
can  be  written  as: 

\au\  =  a+u*  -a~u~  =  A  (25) 


and 

[Pu,}~iru;-p-u;=B  (26) 

where  the  superscripts,  “+”  and  represent  the  variables  and  constants  at  the  right  and  left  sides 
of  the  interface  T ,  respectively.  The  constants,  a+ ,  a~ ,  p* ,  (T ,  A ,  and  B ,  are  known 
constants  determined  by  the  nature  of  the  equation  being  computed.  In  an  actual  two-phase  flow 
problem,  the  jump  conditions  can  be  formulated  such  that  a* ,  a~ ,  /T ,  /T  are  dimensionless 
constants. 


A  grid  point  is  called  a  regular  point  if  the  finite  difference  formulas  at  this  point  only  involve  grid 
points  on  the  same  side  of  the  interface.  Otherwise,  it  is  an  irregular  point.  If  grid  point  i  is  a 
regular  point  (without  the  interface),  finite  difference  approximation  of  an  arbitrary  order  can  be 
easily  derived  by  a  Taylor  series  expansion  or  by  a  polynomial  interpolation.  For  example,  the 

second  and  fourth  order  central  difference  approximations  to  ( d2u/dx2)i  are: 


and 


ydx2l 


2m, +M, 


7+1 


+  0(h 2) 


(27) 


f d2u ' 


v t&2 ), 


_  -u,_2  + 1 6u, -  30m,  + 1  6um  -  ui+2  ^  4 

12  h2 


(28) 
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Therefore,  in  the  case  of  0(h2)  approximation  in  the  regular  points,  /  and  /  +  1  (Fig.  19)  are 

irregular  points.  For  the  0(h 4)  approximation,  there  are  four  irregular  points,  from  /-I  to 

i  +  2.  Alternatively,  we  can  treat  only  points  i  and  /  +  1  as  irregular  points,  while  4th  order 
one-sided  approximation  is  used  for  points  i  — 1  and  i  +  2.  The  latter  approach  is  used  in  this 
paper. 

For  irregular  point  i  shown  in  Fig.  19,  the  use  of  Eq.  (27)  will  lead  to  large  computational  errors 

because  of  the  discontinuity  of  u  at  T .  LeVeque  and  Li  [22]  derived  a  locally  0(h) 

approximation  at  i  by  adding  a  correction  term  to  the  three-point  stencil.  The  difference  formula 
and  the  correction  term  are  then  determined  by  a  Taylor  expansion  at  the  interface.  In  order  to  reach 

0(h)  approximation,  it  is  necessary  to  know  the  jump  conditions  involving  u ,  ux  and  ua .  The 

first  two  jump  conditions  can  be  obtained  easily.  The  jump  conditions  for  the  second  derivatives  can 
be  obtained  by  taking  derivatives  of  the  lower  order  jump  conditions  together  with  the  differential 
equation. 

Therefore,  the  high-order  immersed  interface  method  presented  below  achieves  a  high-order 
approximation  at  the  irregular  point  i  by  imposing  the  two  jump  conditions  given  by  Eqs.  (25)  and 
(26)  only.  Instead  of  using  more  and  more  jump  conditions  of  higher  order  derivatives  to  achieve 
higher  order  accuracy  at  the  irregular  point  as  done  in  the  original  IIM  method,  we  use  more  and 
more  grid  points  on  both  sides  of  the  interface  so  that  arbitrary  order  approximation  can  be  achieved 
while  only  the  two  jump  conditions  (25)  and  (26)  are  imposed.  General  difference  formulas  for 

(dul  dx)j  and  ( d2ul  dx2)t  in  explicit  form  can  be  derived  by  a  matched  polynomial  expansion,  so 

that  they  are  problem  independent. 


Finite  difference  approximation  for  ( du/dx)i  and  ( d2u/dx2)j  at  the  irregular  point  /  is 

considered  by  using  a  stencil  of  n  points  on  the  left  of  T  and  m  points  on  the  right  (Fig.  2).  The 
order  of  the  approximation  increases  with  the  increasing  values  of  n  and  m  .  In  order  to  have  a 
uniform  accuracy,  it  is  desirable  to  have  the  same  number  of  points  on  both  sides  of  the  interface, 
i.e. 

n  =  m  (29) 

Since  n  and  m  can  be  different  in  some  special  circumstances,  we  derive  the  general  formulas 
by  assuming  arbitrary  values  of  n  and  m  .  The  finite  difference  formulas  for  i  can  be  derived 

by  a  Taylor  series  expansion  with  a  as  a  parameter.  The  case  of  locally  0(h)  approximation  for 
(d2u/ dx1)t  at  interface  is  considered  first  as  an  example  below.  A  general  formula  for  arbitrary 
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values  of  n  and  m  is  derived  afterward. 


6.3  Test  and  applications 

The  new  method  can  be  of  arbitrarily  high-order  accuracy  and  it  is  simple  to  be  applied  to  practical 
two-phase  flow  problems  by  requiring  only  the  physical  jump  conditions  for  variables  and  first 
derivatives.  It  also  has  the  advantage  that  the  finite  difference  formulas  at  irregular  points  are 
expressed  in  an  explicit  form  so  that  they  can  be  applied  to  difference  problems  without 
modifications.  Six  versions  of  the  new  method  of  up  to  fourth  order  accuracy  have  been  tested  for 
both  one  and  two-dimensional  model  equations.  The  numerical  results  show  that  they  can  produce 
very  accurate  results  for  elliptic  equations  with  embedded  interfaces.  Details  of  test  results  have 
been  published  in  Publications:  [2],  [12]. 
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